Escaped Muons

Plot list (for 3D plots scroll down from the legend; you can rotate and zoom the 3D plots with your mouse):
Volume Category Histograms | Where muons die (coffee ring plots) | Muon ranging out in iron yoke | Muon scatters in inflector | Distance distribution in Ring iron | Energy loss in volumes | Number of steps in volumes | Another muon in the iron | Stored muon | Coherent Betatron Oscillation

1 Introduction

The goal of this document is to examine “escaped muons” that do not decay in the storage region and do not enter a calorimeter. Some of these muons may leave the storage ring altogether.

I generated art files of simulated events using mdc2 with different magnetic fields turned on or off (scenarios). Note: These files are still not the latest version of MDC2. Nathan is tuning up the injection gun. So right now the storage efficiency I see (0.5%) is still low. Should get better with the next round. We want to examine,

• Muons that are eaten by the iron yoke and muons that escape the world volume entirely.
• Where do muons escape the world (3D and heat map plots?)
• Comparison of the different scenarios

# Load necessary libraries
library(reticulate)  # Access to python
library(stringr)   # string functions
library(parallel)  # Parallel processing (built-in to R)
library(dplyr)     # data analysis
library(purrr)     # Functional programming
library(readr)     # Read formats
library(ggplot2)   # ploting
library(rgl)       # 3D plots
library(gridExtra) # Arranging plots
library(testthat)  # testing
#
library(readGallery) # Read art Gallery files

# Put all readGallery::useDataProduct calls here
useDataProduct('std::vector<gm2truth::TrackingActionArtRecord>')
useDataProduct('std::vector<gm2ringsim::GeantTrackRecord>')
useDataProduct('std::vector<gm2truth::GhostDetectorArtRecord>')

# To restore most of the environment
trackingActionDF <- read_rds('trackingActionDF.rds')
trackingActionNDF <- read_rds('trackingActionNDF.rds')
volNameDF <- read_rds('volNameDF.rds')
ghCylDF <- read_rds('ghCylDF.rds')
ghWldDF <- read_rds('ghWldDF.rds')

2 Load the samples

Samples are located on the Fermilab dCache in the directory /pnfs/GM2/scratch/users/lyon/arr_20170321.

# Function to properly alter paths to accomodate XRootD
#    e.g. convert /pnfs/BLA --> root://fndca1.fnal.gov/pnfs/fnal.gov/usr/BLA
xrootd_ify <- function(aPath) paste0('root://fndca1.fnal.gov/pnfs/fnal.gov/usr', str_replace(aPath, '/pnfs', ''))

Determine the locations of the data files

# Determine locations of our files
system("ifdh ls '/pnfs/GM2/scratch/users/lyon/arr_20170321/*/*.root' | grep ARR", intern=T) %>% sort() -> arrFiles

arrFiles %>% xrootify() -> arrXrdFiles
arrXrdFiles

 [1] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root"   
 [2] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root"     
 [3] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root"  
 [4] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root"
 [5] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root"     
 [6] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root"      
 [7] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root"   
 [8] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root"     
 [9] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root"  
[10] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root"
[11] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root"     
[12] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root"      

write_rds(arrFiles, 'arrFiles.rds')

Pull out the scenaios (these will be in the same order as the file names)

scenarios <- str_match(arrFiles, 'arr_20170321/.+unified_(.+)_cyl')[,2]
scenarios %>% unique()

[1] "everything"    "noDipole"      "noInflector"   "noKickerQuads" "noKicker"      "noQuads"      

Let’s look at the contents of these files (they’re all the same, so we’ll just do one) gm2 v7_04_00 will have product_sizes_dumper.

stringToRun <- str_interp(
  "ssh lyon@gm2gpvm04.fnal.gov 'source /cvmfs/gm2.opensciencegrid.org/prod7/g-2/setup
   setup gm2 v7_04_00 -q prof ; 
   product_sizes_dumper ${aFile} | grep ::'", 
  list(aFile = arrFiles[1] %>% str_replace('fnal.gov/usr/',''))
)
#
system(stringToRun)

        1410762186      0.456 gm2truth::MagnetIronAndCryoArtRecords_artg4_magnetIronAndCryostat_mdc2arr.
         660570363      0.214 gm2truth::TrackingActionArtRecords_artg4__mdc2arr.
         637307490      0.206 gm2truth::TrajectoryArtRecords_artg4__mdc2arr.
         240738334      0.078 gm2ringsim::GeantTrackRecords_artg4_KeepStepsAction_mdc2arr.
         118035827      0.038 gm2truth::RingArtRecords_artg4_Ring_mdc2arr.
          12078805      0.004 gm2truth::GhostDetectorArtRecords_artg4_GhostCylinderDetector_mdc2arr.
           4482185      0.001 art::RNGsnapshots_randomsaver__mdc2arr.
           4086898      0.001 gm2truth::RingTrackingPlaneArtRecords_artg4_RingTrackingPlanes_mdc2arr.
           1373178      0.000 gm2truth::GhostDetectorArtRecords_artg4_GhostNearWorldDetector_mdc2arr.
           1185772      0.000 gm2truth::PhaseSpaceArtRecord_artg4__mdc2arr.
             36971      0.000 art::TriggerResults_TriggerResults__mdc2arr.

3 Escaping muons from the Tracking Action data

The tracking action data has birth and death data for every Geant track.

3.1 Load Tracking Action data

We’ll capture this information for the primary muon and its first generation daughters. See gm2dataproducts/mc/actions/track/TrackingActionArtRecord.hh.

We need a readGallery reader python class. Generate with,

readerClassSkel('gm2truth::TrackingActionArtRecord', writeFile = 'trackingActionReader.py')

Here is the reader,

read_file('trackingActionReader.py') %>% cat

from readGallery import GalleryReaderBase  # Necessary for the base class

class TrackingActionArtRecordReader(GalleryReaderBase):
  def __init__(self, inputTag):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID', 
                  'parentTrackID', 'turn', 'volumeUID', 'status',  'p', 'e', 'x', 'y', 'z', 'px', 'py', 'pz']

  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))

  def fill(self, ROOT, ev):

    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for 
                                                      # gm2truth::TrackingActionArtRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2truth::TrackingActionArtRecord
      for e in p:                             # Loop over elements and fill
        
        # Only accept the primary muon death and its first generation daughters birth
        if (e.trackID == 1 and e.status == 1) or (e.parentTrackID == 1 and e.status == 0):
        
          self.vals.append(
            [ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID, 
              e.parentTrackID, e.turn, e.volumeUID, e.status, e.p, e.e, 
              e.x, e.y, e.z, e.px, e.py, e.pz ])

    return True

Make the reader class and object

trackingActionReaderC <- createReaderClass_from_file('trackingActionReader.py')$TrackingActionArtRecordReader

We have many files to process. Let’s do this reading in parallel.

We need to create a reader object per parallel job

trackingActionReaders <- map(1:length(arrXrdFiles), function(i) trackingActionReaderC(artInputTag("artg4")))

Load the data

# For some reason these jobs seem to run serially
jobs <- lapply(1:length(arrXrdFiles), function(i) getGalleryData(arrXrdFiles[i], trackingActionReaders[[i]]) %>%
                                                    mcparallel(name=i))
trackingActionReturns <- mccollect(jobs)

Let’s merge all of the output.

trackingActionDF <- map_df(1:length(trackingActionReaders), 
                           function(i) galleryReader_df(trackingActionReaders[[i]]) %>% 
                             mutate(scenario=scenarios[i],fileEntry=i))

nrow(trackingActionDF)

[1] 488953

trackingActionDF %>% head(100)  # Let's not dump the entire data frame into the HTML page

# Write this out so we don't need to do the above again
write_rds(trackingActionDF, 'trackingActionDF.rds')

How many muons per file?

trackingActionDF %>% filter(trackID==1) %>% group_by(scenario, fileEntry) %>% tally()

3.1.1 Add Volume Names

We need to turn the volume IDs into volume names. The volume IDs change for each file. We can run a small art job to determine the volume ID to name tables.

runForVolIDString <- function(i) {
  str_interp(
    "PVS_CSVOUT=${csvout}_${i}_volNames.csv gm2 -c ${fclPath}/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 ${aFile}",
    list( csvout=scenarios[i], i=i, fclPath=Sys.getenv("MRB_BUILDDIR"), aFile=arrXrdFiles[i]) 
  )
}
runForVolIDStrings <- sapply(1:length(arrXrdFiles), runForVolIDString)
runForVolIDStrings

 [1] "PVS_CSVOUT=everything_1_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root"       
 [2] "PVS_CSVOUT=noDipole_2_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root"           
 [3] "PVS_CSVOUT=noInflector_3_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root"     
 [4] "PVS_CSVOUT=noKickerQuads_4_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root" 
 [5] "PVS_CSVOUT=noKicker_5_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root"           
 [6] "PVS_CSVOUT=noQuads_6_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root"             
 [7] "PVS_CSVOUT=everything_7_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root"       
 [8] "PVS_CSVOUT=noDipole_8_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root"           
 [9] "PVS_CSVOUT=noInflector_9_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root"     
[10] "PVS_CSVOUT=noKickerQuads_10_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root"
[11] "PVS_CSVOUT=noKicker_11_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root"          
[12] "PVS_CSVOUT=noQuads_12_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root"            

Let’s run the art jobs in parallel

jobs <- lapply(1:length(arrXrdFiles), function(i) mcparallel(system( runForVolIDStrings[i], intern=T ), name=i))
res <- mccollect(jobs)

and load the CSV files.

jobs <- lapply(1:length(arrFiles), function(i) 
  str_interp("${scn}_${i}_volNames.csv", list(scn=scenarios[i], i=i)) %>%
    read_csv(col_names=c("volumeUID", "volName")) %>% mcparallel(name=i) )
#
volNameTables <- mccollect(jobs)
volNameDF <- map_df(1:length(volNameTables), function(i) 
                      volNameTables[[i]] %>% mutate(scenario=scenarios[i],
                                                    fileEntry=i))
nrow(volNameDF)

[1] 1014

volNameDF %>% head(100)

write_rds(volNameDF, 'volNameDF.rds')

Now we merge with the tracking action data frame

trackingActionDF %>% inner_join(volNameDF) %>% 
  select(scenario, fileEntry, eventEntry, volName, everything()) -> trackingActionNDF

Joining, by = c("fileEntry", "volumeUID", "scenario")

nrow(trackingActionNDF)

[1] 488953

trackingActionNDF %>% head(100)

3.2 How many muons decay in the storage region?

These muons decay in the storage region (not quite the same as “stored”),

trackingActionNDF %>% filter(trackID==1, volName=='ArcSection[00]') %>% 
  group_by(scenario, fileEntry) %>% tally()

Interesting! So, from this, only 0.575% decay in the storage region for the “all on” (or everything) sample. That seems rather low. But continuing.

write_rds(trackingActionNDF, 'trackingActionNDF.rds')

3.3 Where do muons go to die?

trackingActionNDF %>% filter(trackID == 1) %>% group_by(scenario, volName) %>% tally()

It may be nice to categorize the volume names by area (e.g. “Arc”, “bellows”, “Ring”, “Inflector”).

# Pull out the first word in camel case (see http://stackoverflow.com/questions/29916065/how-to-do-camelcase-split-in-python)
camelCaseSplit <- function(s) str_split(s, '(?<=[a-z])(?=[A-Z])|(?<=[A-Z])(?=[A-Z][a-z])')

# Test the above
expect_equal(camelCaseSplit("RingPoleTipUpper")[[1]], c("Ring", "Pole", "Tip", "Upper"))
expect_equal(camelCaseSplit("ringPoleTipUpper")[[1]], c("ring", "Pole", "Tip", "Upper"))
expect_equal(camelCaseSplit("world")[[1]], c("world"))
expect_equal(camelCaseSplit("ringPole")[[1]], c("ring", "Pole"))

# We really just want the first word
takeFirstCamel <- function(s) camelCaseSplit(s) %>% map_chr(function(s) s[1])
#
# Test
takeFirstCamel("RingPoleTipUpper") %>% expect_equal("Ring")
takeFirstCamel("ringPoleTipUpper") %>% expect_equal("ring")
takeFirstCamel("world") %>% expect_equal("world")
takeFirstCamel("ringPole") %>% expect_equal("ring")
takeFirstCamel(c("ringPoleTipUpper", "world")) %>% expect_equal(c("ring", "world"))

Don’t confuse RingTrackingPlane with the iron (also starts with Ring)

makeVolCategory <- function(s) ifelse( str_detect(s, 'RingTracking'),
                                       "RngTrkPlane",
                                       takeFirstCamel(s))

makeVolCategory(c("RingBottom", "RingTrackingPlane[8]", "world")) %>% expect_equal(c("Ring", "RngTrkPlane", "world"))

trackingActionNDF %>% mutate(volCategory = volName %>% makeVolCategory()) -> trackingActionNVDF
trackingActionNVDF %>% head(100)

Show count by category

trackingActionNVDF %>% filter(trackID == 1) %>% group_by(scenario, volCategory) %>% tally() -> volCatTally
volCatTally %>% mutate(perc = n/20000 * 100) -> volCatTally
volCatTally

Let’s plot!

trackingActionNVDF %>% filter(trackID ==1) %>% ggplot(aes(x=volCategory, group=scenario)) + 
  geom_bar()

ggplot(volCatTally, aes(x = volCategory, y=perc, fill=scenario)) + 
  geom_bar(stat="Identity", width=0.5, position="dodge")

Can’t really see non-Ring/world ones. Let’s split them out.

# Function to make this easy
plotit <- function(d) ggplot(d, aes(x = volCategory, y=perc, fill=scenario) ) + 
  geom_bar(stat="Identity", width=0.5, position="dodge") + 
  xlab("Volume category") + ylab("Percent")
volCatTally %>% filter(volCategory %in% c("Ring", "world")) %>% plotit -> p1
volCatTally %>% filter(! volCategory %in% c("Ring", "world")) %>% plotit -> p2
grid.arrange(p1, p2)

Let’s see where they die

# See https://rpubs.com/hadley/97970 for how to wrap a multipart ggplot2 component
plotCommon <- function () {
  list(
    facet_wrap(~scenario),
    guides(col = guide_legend(override.aes = list(size=5, alpha=1))),
    scale_color_discrete(name="Volume\nCategory"),
    labs(x="z (mm)", y="x (mm)"),
    theme_minimal()
  )
}

The “Coffee Ring” plots

trackingActionNVDF %>% filter(trackID == 1) %>% 
  ggplot( aes(x=z, y=x, color=volCategory)) + 
    geom_point(alpha=0.1) + 
    ggtitle('Where Muons Die') +
    plotCommon()

Here it is without the ring losses, which dominate

trackingActionNVDF %>% filter(trackID == 1, volCategory != 'Ring') %>% 
  ggplot( aes(x=z, y=x, color=volCategory)) + 
    geom_point(alpha=0.5) + 
    ggtitle("Where Muons Die", subtitle = "Ring escapees excluded for clarity") +
    plotCommon()

3.4 Does it make sense for so many muons to stop in the iron?

We have 3 GeV muons. Does it make sense that the vast majority of them stop in the iron? A table shows the Continuous Slow Down Approximation range of muons (CSDA). For a 3 GeV muon in iron, the CSDA is \(1.825 \times 10^{3}\) g/cm\(^2\). The density \(\rho\) is 7.874 g/cm\(^3\). So, where \(R\) is range,

\[R = \text{CSDA}/\rho\] Range is thus 2.32 m. This is bigger than the width of the iron, but remember that the muons come in at an oblique angle. We’ll have to prove this.

4 How far do muons go in the iron?

We have the Geant stepping information, so we should be able to figure out how far the muons go in the iron.

Here is the Gallery reader,

read_file('GeantTrackReader.py') %>% cat

from readGallery import GalleryReaderBase  # Necessary for the base class

class GeantStepReader(GalleryReaderBase):
  def __init__(self, inputTag, eventsToFill = []):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'trackID', 'globalStepNum', 'totalEnergyDeposit', 
                  'stepLength', 'volumeUID', 
                  'x', 'y', 'z', 'px', 'py', 'pz']
    self.eventsToFill = eventsToFill
    print self.eventsToFill

  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2ringsim.GeantTrackRecord))

  def fill(self, ROOT, ev):

    # Do we care about this event - if not then don't bother looking?
    if self.eventsToFill and not ev.eventEntry() in self.eventsToFill:
      return True
      
    print 'Reading entry %s' % ev.eventEntry()

    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for gm2ringsim::GeantTrackRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2ringsim::GeantStep
      for e in p:                             # Loop over elements and fill
      
        # Get the steps
        for f in e.steps():

          self.vals.append(
            [ ev.fileEntry(), ev.eventEntry(), e.trackID(), f.globalStepNum(), 
            f.totalEnergyDeposit(), f.stepLength(), f.volumeUID(), 
            f.pos()[0], f.pos()[1], f.pos()[2], 
            f.p()[0], f.p()[1], f.p()[2] ])
            
        # If we're on the last event, don't bother reading any more
        if self.eventsToFill and ev.eventEntry() == self.eventsToFill[-1]:
          return False

    return True

geantStepReaderC <- createReaderClass_from_file('GeantTrackReader.py')$GeantStepReader

4.1 Muon distance in iron for one event

Let’s find an event where the muon did in the iron.

trackingActionNVDF %>% filter(trackID == 1, volName == 'RingYokeBottom') %>% head(100)

Event 4 looks good.

geantStepOne <- geantStepReaderC(artInputTag('artg4:KeepStepsAction'), tuple(as.integer(4)))

(4,)

getGalleryData(arrXrdFiles[1], geantStepOne)

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Reading entry 4
Timings: Overall time = 5.234 s
Time per event: min=0.000   mean=0.206   max=1.027

gsOneDF <- galleryReader_df(geantStepOne)

Merge in volume names

gsOneDF %>% mutate(fileEntry = 1) %>%
  inner_join(volNameDF) %>% 
  select(volName, everything()) -> gsOneDF

Joining, by = c("fileEntry", "volumeUID")

gsOneDF

Let’s make a 3D plot of this path

ring <- cylinder3d( center=rbind(c(0,-90, 0), c(0,90,0)),
                    radius=7112,
                    sides=20, closed = F)

Let’s do volume categories again.

gsOneDF %>% mutate(volCategory = volName %>% makeVolCategory %>% as.factor) -> gsOneDF

Let’s plot the path. Note that I have to make the legend separately since legend3d looks awful.

# lengend3d looks terrible - so do a regular legend
plot(1, type='n', axes=FALSE, ann=FALSE)
with(gsOneDF,
  legend("top", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)

clear3d()
with(gsOneDF,{
     plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

How much energy loss and distance traveled in each volume category

gsOneDF %>% group_by(volCategory) %>% summarize(totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n())

Ok - I believe this (note that the y scale is very small compared to the x,z scales - so the muon acutally travels quite far in x & z). The muon travels about 2.4m in the iron and ranges out. Physics (and Geant) work!

Let’s try another one

geantStepOneA <- geantStepReaderC(artInputTag('artg4:KeepStepsAction'), tuple(as.integer(10)))

(10,)

getGalleryData(arrXrdFiles[1], geantStepOneA)

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Reading entry 10
Timings: Overall time = 2.566 s
Time per event: min=0.000   mean=0.090   max=0.993

gsOneADF <- galleryReader_df(geantStepOneA)
gsOneADF %>% mutate(fileEntry = 1) %>%
  inner_join(volNameDF) %>% 
  mutate(volCategory = volName %>% makeVolCategory %>% as.factor)%>% 
  select(volName, volCategory, everything()) -> gsOneADF

Joining, by = c("fileEntry", "volumeUID")

gsOneADF

gsOneADF %>% group_by(volCategory) %>% summarize(totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n()) -> gsOneADFSum
gsOneADFSum

Check sums

gsOneADFSum %>% summarize(sum(totalELoss_MeV))

gsOneADF %>% summarize(sum(totalEnergyDeposit))

Sanity check passes!

# lengend3d looks terrible - so do a regular legend
plot(1, type='n', axes=FALSE, ann=FALSE)
with(gsOneADF,
  legend("top", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)

clear3d()
with(gsOneADF,{
     plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
            zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

This one scatters in the inflector.

4.2 Muon distance for many events

Can we show how far muons go in various materials?

Let’s collect data for 100 events.

We have a reader that can do n events.

read_file('GeantTrackReaderN.py') %>% cat

from readGallery import GalleryReaderBase  # Necessary for the base class

class GeantStepReaderN(GalleryReaderBase):
  def __init__(self, inputTag, startEvent=100, nEvents=10):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'trackID', 'globalStepNum', 'globalTime', 'totalEnergyDeposit', 
                  'stepLength', 'volumeUID', 
                  'x', 'y', 'z', 'px', 'py', 'pz']
    self.startEvent = startEvent
    self.endEvent   = startEvent + nEvents
    
  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2ringsim.GeantTrackRecord))

  def fill(self, ROOT, ev):

    # Do we care about this event - if not then don't bother looking
    if ev.eventEntry() < self.startEvent:
      return True
      
    if ev.eventEntry() >= self.endEvent:
      return False

    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for gm2ringsim::GeantTrackRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2ringsim::GeantStep
      for e in p:                             # Loop over elements and fill
      
        # Get the steps
        for f in e.steps():

          self.vals.append(
            [ ev.fileEntry(), ev.eventEntry(), e.trackID(), f.globalStepNum(), f.globalTime(), 
            f.totalEnergyDeposit(), f.stepLength(), f.volumeUID(), 
            f.pos()[0], f.pos()[1], f.pos()[2], 
            f.p()[0], f.p()[1], f.p()[2] ])

    return True

geantStepReaderNC <- createReaderClass_from_file('GeantTrackReaderN.py')$GeantStepReaderN

geantStepN <- geantStepReaderNC(artInputTag('artg4:KeepStepsAction'), 1000, 1000)

getGalleryData(arrXrdFiles[1], geantStepN)

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Timings: Overall time = 51.662 s
Time per event: min=0.000   mean=0.024   max=0.452

gsNDF <- galleryReader_df(geantStepN)
gsNDF %>% mutate(fileEntry = 1,
                 p = sqrt(px*px + py*py + pz*pz)) %>%
  inner_join(volNameDF) %>% 
  mutate(volCategory = volName %>% makeVolCategory %>% as.factor)%>% 
  select(volName, volCategory, p, everything()) -> gsNDF

Joining, by = c("fileEntry", "volumeUID")

nrow(gsNDF)

[1] 297501

gsNDF %>% head(100)

Let’s do a distribution of distance in the different volumes

gsNDF %>% group_by(eventEntry, volCategory) %>% summarize(
                                                totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n()) -> gsNDFSum
gsNDFSum

gsNDFSum %>% filter(volCategory == 'Ring') %>% 
  ggplot(aes(x=totalDist_meters)) + geom_histogram(bins=50) + 
    labs(x = 'Total distance in Ring material (m)', y='count', title='Muon in Ring Material')

Let’s look at this across the volume categories for energy loss

gsNDFSum %>% ggplot(aes(x=totalELoss_MeV)) + geom_histogram(bins=50) +
  facet_wrap(~volCategory, scales = "free") + 
  labs(x = 'Total energy loss (MeV)', y='count', title='Energy loss in material',
       subtitle = 'Note the different scales')

Let’s get an energy loss curve

gsNDF %>% ggplot(aes(x=p, y=totalEnergyDeposit)) + geom_point() + facet_wrap(~volCategory) + 
  labs(x="Muon Momentum (MeV/c)", y="Total Energy loss (MeV))")

Well, not sure what all this means. Let’s look for a particular event…

gsNDF %>% filter(eventEntry == 1300) %>% ggplot(aes(x=p, y=totalEnergyDeposit)) + geom_point() + facet_wrap(~volCategory) + 
  labs(x="Muon Momentum (MeV/c)", y="Total Energy loss (MeV))", title='Energy loss vs. Momentum for muons in event 1300' )

How many steps do we take?

gsNDFSum %>% ggplot(aes(x=nSteps)) + geom_histogram(bins=50) +
  facet_wrap(~volCategory, scales = "free") + 
  labs(x = 'Number of steps', y='count', title='Number of steps in materials',
       subtitle = 'Note the different scales')

Let’s look at two events

4.2.1 Muon with long distance in iron

gsNDFSum %>% filter(volCategory == 'Ring', totalDist_meters > 3)

Let’s look at event 1758

gsNDF %>% filter(eventEntry == 1758) -> gsNDF1758
gsNDF1758

gsNDF1758 %>% group_by(eventEntry, volCategory) %>% summarize(
                                                totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n()) -> gsNDF1758Sum
gsNDF1758Sum

Let’s plot it!

gsNDF1758 %>% mutate(volCategory = as.factor(as.character(volCategory))) -> gsNDF1758

# lengend3d looks terrible - so do a regular legend
plot(0, type='n', axes=FALSE, ann=FALSE)
with(gsNDF1758,
  legend("bottom", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)

clear3d()
with(gsNDF1758,{
     plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
            zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

Looks like this one started out in a strange place.

4.2.2 Muon with many turns

gsNDFSum %>% filter(volCategory == 'Arc', totalDist_meters > 870)

Let’s try event 1509

gsNDF %>% filter(eventEntry == 1509) -> gsNDF1509
gsNDF1509

gsNDF1509 %>% group_by(eventEntry, volCategory) %>% summarize(
                                                totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n()) -> gsNDF1509Sum
gsNDF1509Sum

Let’s plot it!

gsNDF1509 %>% mutate(volCategory = as.factor(as.character(volCategory))) -> gsNDF1509

# lengend3d looks terrible - so do a regular legend
plot(0, type='n', axes=FALSE, ann=FALSE)
with(gsNDF1509,
  legend("bottom", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)

clear3d()
with(gsNDF1509,{
     plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
            zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

Looks like a nice stored muon!! Can we see CBO?

gsNDF1509 %>% mutate(r = sqrt(x*x+z*z)) -> gsNDF1509

p1 <- gsNDF1509 %>% filter(volCategory %in% c('RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=r)) + geom_line() +
  labs(x='Global time (ns)', y='Radius (mm)', title='Radial CBO', subtitle="Measured by Ring Tracking Planes")
#
p2 <- gsNDF1509 %>% filter(volCategory %in% c('Arc','RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=r)) + 
  geom_line() + 
  labs(x='Global time (ns)', y='Radius (mm)', subtitle="Measured by Arc hits and Ring Tracking Planes")
#
grid.arrange(p1, p2)

We can even see the kick!

Here we’ll plot the points too.

p1 <- gsNDF1509 %>% filter(volCategory %in% c('RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=y)) + geom_line() +
  geom_point() +
  labs(x='Global time (ns)', y='y (mm)', title='Vertical CBO', subtitle="Measured by Ring Tracking Planes")
#
p2 <- gsNDF1509 %>% filter(volCategory %in% c('Arc','RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=y)) + 
  geom_line() + geom_point() + 
  labs(x='Global time (ns)', y='y (mm)', subtitle="Measured by Arc hits and Ring Tracking Planes")
#
grid.arrange(p1, p2)

5 “Ghost” cylinder at vacuum chamber

We want to examine muon escapees that leave the vacuum chamber. To measure this feature, we have a ghost cylinder at the outer wall of the vacuum chamber as well as above and below the storage region (before the muon would hit iron). That should catch every muon that doesn’t get stored and makes it out of the storage region. The data is gm2truth::GhostDetectorArtRecords_artg4_GhostCylinderDetector.

There is also a ghost cylinder just on the inside of the world cube to see muons that would actually leave the ring (and make it through the iron). The data is gm2truth::GhostDetectorArtRecords_artg4_GhostNearWorldDetector.

5.1 Load the Ghost cylinder data

See gm2dataproducts/mc/ghostdetectors/GhostDetectorArtRecord.hh .

readerClassSkel('gm2truth::GhostDetectorArtRecord', writeFile = 'GhostDetectorReader.py')

read_file('GhostDetectorReader.py') %>% cat

from readGallery import GalleryReaderBase  # Necessary for the base class

class GhostDetectorArtRecordReader(GalleryReaderBase):
  def __init__(self, inputTag):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'particleID', 
                  'trackID', 'parentTrackID', 
                  'x', 'y', 'z', 
                  'px', 'py', 'pz', 
                  'sx', 'sy', 'sz', 'time', 'volumeName']

  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.GhostDetectorArtRecord))

  def fill(self, ROOT, ev):

    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for gm2truth::GhostDetectorArtRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2truth::GhostDetectorArtRecord
      
      for e in p:                             # Loop over elements and fill

        self.vals.append(
          [ ev.fileEntry(), ev.eventEntry(), e.particleID, e.trackID, e.parentTrackID, 
            e.position.x(), e.position.y(), e.position.z(), 
            e.momentum.x(), e.momentum.y(), e.momentum.z(), 
            e.spin.x(), e.spin.y(), e.spin.z(), e.time, e.volumeName ])

    return True

ghostDetectorReaderC <- createReaderClass_from_file('GhostDetectorReader.py')$GhostDetectorArtRecordReader

Let’s read both ghost cylinders

ghCylReader <- ghostDetectorReaderC(artInputTag('artg4:GhostCylinderDetector'))
ghWldReader <- ghostDetectorReaderC(artInputTag('artg4:GhostNearWorldDetector'))

Load for one file

getGalleryData(arrXrdFiles, c(ghCylReader, ghWldReader))

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Closing file, read 16740159 bytes in 100 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root
Closing file, read 11886702 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root
Closing file, read 14680905 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root
Closing file, read 16282559 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root
Closing file, read 16397132 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root
Closing file, read 16580366 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root
Closing file, read 16737803 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root
Closing file, read 11850233 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root
Closing file, read 14723528 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root
Closing file, read 16377209 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root
Closing file, read 16397059 bytes in 66 transactions
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root
Closing file, read 16814838 bytes in 66 transactions
Timings: Overall time = 252.037 s
Time per event: min=0.000   mean=0.000   max=7.044

Load the data into R

ghCylDF <- galleryReader_df(ghCylReader)
ghWldDF <- galleryReader_df(ghWldReader)

Look at what we got

nrow(ghCylDF) %>% format(big.mark=",")

[1] "2,374,697"

nrow(ghWldDF) %>% format(big.mark=",")

[1] "197,905"

ghCylDF %>% head(100)

ghWldDF %>% head(100)

write_rds(ghCylDF, 'ghCylDF.rds')  # Write out just in case, since these are expensive to extract
write_rds(ghWldDF, 'ghWldDF.rds')

Fill in the scenario and pdg

pdgs <- c('e-'= 11, 'nu_e'= 12, 'mu-'= 13, 'nu_mu'= 14, 'tau-'= 15, 'nu_tau'= 16,
          'e+'=-11, 'anti_nu_e'=-12, 'mu+'=-13, 'anti_nu_mu'=-14, 'tau+'=-15, 'anti_nu_tau'=-16,
          'gam'=22, 'p'=2212, 'n'=2112, 'anti-p'=-2212, 'anti-n'=-2112)
pdgs

         e-        nu_e         mu-       nu_mu        tau-      nu_tau          e+   anti_nu_e         mu+  anti_nu_mu 
         11          12          13          14          15          16         -11         -12         -13         -14 
       tau+ anti_nu_tau         gam           p           n      anti-p      anti-n 
        -15         -16          22        2212        2112       -2212       -2112 

# We need to be a little careful, since the scenario names are repeated (we have two files per scenario). 
# We can't make factor levels out of that. 
# The file entries with 0-5 are easy. 6-11 are harder
addScenarioFactor <- function(f) factor(if_else(f >= 6, as.integer(f-6), f), 
                                        levels=0:(length(unique(scenarios))-1), labels=unique(scenarios))
#
mutateThis <- function(df) mutate(df, 
                       scenario = addScenarioFactor(fileEntry),
                            pdg = factor(particleID, levels=pdgs, labels=names(pdgs)),
                              r = sqrt(x*x+z*z),
                            phi = atan2(z, x) %>% if_else(. > 0, ., 2*pi+ .) * (360.0/(2*pi)) )
#
ghCylDF %>% mutateThis -> ghCylDF
#
ghWldDF %>% mutateThis -> ghWldDF

5.2 Losses detected by the cylinder ghost detector

The “ghost cylinder” virtual detector in the parallel world sees muon escapees that hit the vacuum chamber outer wall. How many hits per muon do we see?

Pull out just the muon

ghCylDF %>% filter(trackID == 1) -> ghCylMuDF
nrow(ghCylMuDF) %>% format(big.mark=',')

[1] "132,158"

To do this right, we need to make a “key” out of the fileEntry and eventEntry columns.

ghCylMuDF %>% mutate( flevKey = unite ...  )

ghCylMuDF %>% group_by(fileEntry, eventEntry) %>% tally()

I’ve made 3D plots of all the scenarios (see below). Interesting to see what fraction get kicked up and down (are the ghost hits on the “ceiling or floor” or on the “wall”). We can tell by plotting the y coordinate…

ghCylDF %>% 
  ggplot(aes(y)) + geom_histogram(bins=50) + scale_y_log10() + facet_wrap(~scenario)

5.2.1 3D plots of all scenarios

Let’s practice with just the everything scenario

ghCylMu_everything_DF <- ghCylMuDF %>% filter(scenario == 'everything')
nrow(ghCylMu_everything_DF) %>% format(big.mark=',')

[1] "23,691"

Let’s plot where these are

plotEscapees <- function(i) {
  with( ghCylMuDF %>% filter(scenario == scenarios[i]),
   plot3d(x=x, y=y, z=z, xlim=c(-8000, 8000), ylim=c(-120, 120), zlim=c(-8000, 8000), main=scenarios[i]))
  plot3d(ring, add=T, alpha=0.2)
  view3d(phi=90, theta=-90)
}

clear3d()
plotEscapees(1)
rglwidget()

They’re all lost up or down, not along the side. Let’s try more of the scenarios.

clear3d()
plotEscapees(2)
rglwidget()

No storage ring dipole field: As we expect, the muons don’t go very far.

clear3d()
plotEscapees(3)
rglwidget()

No inflector field: Seemingly lots of muons kicked up and down (need to check that). Many slam into the wall.

clear3d()
plotEscapees(4)
rglwidget()

clear3d()
plotEscapees(5)
rglwidget()

clear3d()
plotEscapees(6)
rglwidget()

Well, this isn’t too illuminating.